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Abstract. We examine the onset of classical topological order in a nearest- neighbor 
kagome ice model. Using Monte Carlo simulations, we characterize the topological 
sectors of the groundstate using a non-local cut measure which circumscribes the 
toroidal geometry of the simulation cell. We demonstrate that simulations which 
employ global loop updates that are allowed to wind around the periodic boundaries 
cause the topological sector to fluctuate, while restricted local loop updates freeze the 
simulation into one topological sector. The freezing into one topological sector can 
also be observed in the susceptibility of the real magnetic spin vectors projected onto 
the kagome plane. The ability of the susceptibility to distinguish between fluctuating 
and non-fluctuating topological sectors should motivate its use as a local probe of 
topological order in a variety of related model and experimental systems. 



1. Introduction 

There has been great interest recently in quantum liquids which fail to find symmetry- 
broken states down to temperatures far below that set by the leading energy scale. 
Rather, such systems occupy an extensive band of low lying states which retain some of 
the fluid-like features of high-temperature phases even in the absence of configurational 
disorder. Magnetic systems have become the canonical systems for study and this 
disordered state is referred to as a quantum spin liquid [1, 2]. Spin liquids, however, 
differ in important ways from their higher-temperature paramagnetic counterparts, due 
to the existence of subtle correlations resulting from highly-constrained collective motion 
of the magnetic moments over the low energy manifold of states. The constraints can 
lead to hidden topological information encoded into the low energy states, and quantum 
spin liquids can be classified by this topological order buried in their disordered quantum 
groundstates [3]. 

The concept of topological order can also be applied to classical systems [4, 5]. A 
classical system with topological order is characterized by a highly degenerate manifold 
of ground states whose configurations are subject to a local topological constraint. 
Imposing the constraint everywhere leads again to information being encoded non- 
locally into each configuration, even in the absence of symmetry breaking, or selection 
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of any particular state down to zero temperature. In this paper, we examine a concrete 
and experimentally realizable example of a topologically constrained system, that of 
"kagome ice" [6, 7]. This phase can be realized to a good approximation in spin 
ice materials such as Dy 2 Ti 2 07 by placing a magnetic field of intermediate strength 
along the [111] crystallographic direction. With this field configuration a magnetization 
plateau occurs that corresponds to a dimensional reduction of the magnetic degrees 
of freedom into the kagome planes perpendicular to the field [8]. Application of the 
ice rules for each tetrahedral unit of spins is equivalent to applying a local topological 
constraint to the perpendicular, in-plane degrees of freedom. Within the approximation 
of nearest neighbour interactions only, the resulting degenerate manifold of states is 
identical to two systems: hard core dimers close-packed on a honeycomb lattice; and, 
the groundstate of an antiferromagnetic Ising spin model on the kagome lattice in the 
presence of a magnetic field. Such discrete systems, subject to a local constraint, 
are characterized by an extensive ground state degeneracy and a gap for excitations 
out of the constrained manifold. These excitations break the constraints leading to 
the proliferation of deconfined topological defects [9]. In this paper, we define and 
review several measures of topological order through a topological number or sector. 
We then use Monte Carlo simulations to study the fluctuations of these quantities 
and the onset of topological order in the kagome ice model. Using local and non- 
local updates, we examine the onset of ergodicity-breaking between different topological 
sectors as one cools into the groundstate phase. We demonstrate that one can discern the 
ergodicity breaking corresponding to the onset of topological order using the magnetic 
susceptibility - a local quantity that is easily measurable and accessible in principle 
in experiments. We show that topological order occurs here uniquely through the 
application of constrained dynamics, which limits collective "loop" moves to clusters of 
microscopic length. We compare our results with a collective paramagnetic system [10] 
without topological constraints, so called "Wills-ice" [11], whose ground state manifold 
maps onto the complete low energy phase space of the Ising antiferromagnet on the 
kagome lattice. 

2. Pseudospin Hamiltonian and Topological Sectors 

We begin with the Ising antiferromagnet on the kagome lattice with Hamiltonian: 

(ij) i 

where we refer to the (j; = ±1 as pseudospin variables and where h is an external field 
in reduced units. The kagome lattice is built of a triangular Bravais lattice with a three 
spin basis which we take to be the spins on a "left" triangle (see Figure 1). In zero 
field the best this highly frustrated system can do in lowering its energy is to assure two 
spins up and one down on each triangle (uud)- or the inverse (ddu). This leads to an 
extensive ground state entropy, S/N = 0.502 ke [12], and although the configurations 
on each triangle are not independent, all six combinations of uud and ddu figure in the 
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ground state manifold. Applying a magnetic field breaks the Z 2 symmetry of uud and 
ddu states but still does not lead to an ordered groundstate. Rather, the groundstate 
entropy is reduced, but remains finite at S/N = 0.108 ke [13]. One can make a direct 
mapping between the antiferromagnet and a 2D ice-like system by defining in-plane spin 
degrees of freedom lying along the axes joining the triangles of the kagome lattice. The 
mapping between the real spins Sj and pseduospins <7j is 

Si = diOi (2) 

where di are the three unit vectors joining a left triangle to a right triangle in Figure 1, 
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which leads to a frustrated ferromagnetic Hamiltonian whose exchange term reads 

(ij) 

Within the real spin description, the reduced manifold of states, with broken Z 2 
symmetry corresponds to either two spins pointing out of and one into each left triangle 
(see the Appendix, Figure 9), or vise versa, while the phase space with full symmetry 
corresponds to a mixture of both. Consider the case with reduced symmetry, with two- 
out and one-in and define an effective field Bj such that B; = 2S« for a spin entering 
a left triangle (and leaving a right triangle) and that B; = Sj for all other spins. The 
constraint imposes that the total flux of B into each triangle is zero; that is, that the 
ground state condition plus broken Z 2 symmetry corresponds to a discrete, divergence 
free condition, V-B = 0. Breaking this symmetry thus leads to an emergent U(l) gauge 




Figure 1. A section of the kagome lattice (blue) corresponds to a 2D plane in the 
lattice of corner sharing tetrahcdra. The white spins are assumed to be static, pinned 
by a [111] external field, and therefore do not represent magnetic degrees of freedom 
in the model Equation (1). Rather, they are intended to motivate the internal field 
term h affecting the kagome spin degrees of freedom. In the spin configuration shown, 
minority spins arc shown in red. 
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field encoded into each disordered ground state configuration [14]. The constrained 
system also maps onto a close-packed dimer system on a honeycomb lattice by placing a 
dimer along the line joining the centres of each pair of triangles connected by a minority 
spin [15]. Equivalent gauge field descriptions have been made in other systems with 
dimers tiling a bipartite lattice [16]. Through the emergent gauge field, excitations 
that break the constraints decompose into pairs of quasi-particles carrying topological 
charge. It has recently been shown, in the related physics of three dimensional (3D) 
spin ice materials, that in addition to the topological properties the charges interact via 
an effective Coulomb law, making them examples of magnetic monopole quasi-particles 
[9]. 

The name "kagome ice" was first given to a magnetic ice model interacting via 
Hamiltonian (1), with access to the complete phase space of low energy states (h = 0) 
by Wills et. al. [11]. We refer to this system as "Wills-ice" to distinguish it from 
the constrained manifold of states with broken Z 2 symmetry, attained by applying a 
magnetic field along the [111] axis of a 3D model system, the nearest neighbour spin 
ice [8, 15]. In this model, the spins sit on the vertices of a pyrochlore lattice of corner 
sharing tetrahedra and lie along axes joining the centres of the tetrahedra, in analogy 
with the two dimensional (2D) spins described above. The field isolates kagome planes 
perpendicular to the [111] plane pinning the central spins of the tetrahedra, as shown in 
Figure 1. The 3D spins are subject to the ice rules with two spins in and two out of each 
unit, leaving the other three spins on each tetrahedra lying on a triangular element of a 
kagome lattice, with projections in the perpendicular plane. The 3D ice rules mean that 
the projected spin components are constrained to be two-in, one-out on each triangle, 
or the inverse, with the choice dictated by the magnetic field direction along ±[111]. 
While the in-plane components map exactly onto our "real spins" in Hamiltonian (4), 
the components parallel to [111] are equivalent to the pseudo spins [15]. We simulate this 
constraint by applying a field h = 2 J to the pseudo spin Hamiltonian, which corresponds 
to the effective field of a static fourth spin with tetrahedral geometry, associated with 
each triangle, h fixes the energy scale for topological defects and the energy scale on 
which these topological ice rules are obeyed. In what follows we shall refer to the 
constrained kagome ice as "[111] ice". 

The consequences of the topological constraints have previously been discussed in 
the context of hard-core dimer coverings on the honeycomb lattice. This order can be 
characterized in several ways, the most well-known being through the use of winding 
numbers, introduced by Rokshar and Kivelson on the frustrated square lattice [17]. 
Their method extends to all bipartite lattices, and is defined via construction of a 
transition graph between any hard-core dimer covering, |C), and a chosen reference state 
\C). Given that the honeycomb lattice is bipartite, it is possible to orient the dimers of 
\C) from the A to B sub-lattice. Orienting the dimers of \C) in the opposite direction 
and then superimposing the two states produces a transition graph consisting of non- 
interacting loops (see Figure 2). On a two-torus the topological sector of \C) can then 
be characterized by the winding numbers w x ,w y , defined as the net number of clockwise 
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loops minus the net number of counterclockwise loops encircling the torus in the x and 
y direction respectively. These winding numbers are bound between —L < w x ,w y < L, 
where L is the linear system size of the underlying kagome lattice. Note that loops of 
microscopic length (smaller than the system size) contribute zero to the winding number, 
as a line circumscribing the boundary cuts all such loops twice, circulating once in each 
direction. 

An alternative way to determine the winding numbers of a particular configuration 
is to avoid constructing the transition graph or defining a reference state - taking 
instead the measurement directly on the spins of the kagome lattice. This can be 
done by considering closed "cuts" along x and y respectively, and counting the number 
of minority elements - dimers within the dimer picture; minority pseudo spins, or real 
spins in the magnetic picture. As we apply a field to the pseudo spins in the "up" 
direction, the minority spins are "down" (real spins are "out" of the "up" triangles in 
Figure 3). The cut number, Hx\c), Hy\c)i is then defined as the number of down pseudo 
spins, or out real spins encountered along the x and y directions, as shown in Figure 3. 
As these spins indicate the position of dimers on the honeycomb lattice, cut numbers can 
be constructed by adding +1 every time a dimer is cut along the path. This definition 
gives only positive numbers, as it corresponds to the Z 2 symmetry being broken in a 
particular way. Placing the field in the opposite direction would isolate the opposite 
half of phase space and would require the definition of negative numbers. Counting 
the number of down pseudo-spins means that the cut numbers remain valid also for 
unconstrained Wills ice; however the down pseudo- spins are not necessarily minority 
and the probability of encountering one along any path approaches 1/2. 

To make this measure equivalent to the winding numbers the constrained [111] 
ice, as previously defined, one must take the cut numbers of the current state and 
subtract the cut numbers of the reference state (Equations 5-6). Note that there 
is a key difference between this treatment on the honeycomb lattice and that on 
the square lattice: on the honeycomb there is no possible reference state with cut 
numbers identically equal to zero, while on the square lattice using columnar state as 



a) b) c) Wy=-2 




Figure 2. Two dimer coverings of the dual honeycomb lattice possible under the ice 
rules, (a) shows the state of interest |C), (b) shows the reference columnar state |C") 
and (c) shows their associated transition graph with winding numbers. 
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Figure 3. Two possible cuts taken along the original kagome lattice with the 
honeycomb dimer lattice superimposed on top. Unfilled circles represent minority 
spins on the original kagome lattice, while the red arrows indicate dimers oriented 
from the A to B sub-lattices within the honeycomb lattice. 



a) b) c) 




Figure 4. Height model covering of the a) state \C), b) the reference state |C") and 
c) of their transition graph. The difference in heights is constant within each loop. 

the reference, Hx\c>), Hy\c) — 0, so Hx\o), Hy\c) = w x ,w y for every possible dimer 
covering. 

w x = Hx\ C ) - Hx\c) (5) 

w y = Hy\ c) - Hy\c) (6) 

A third possibility for the definition of the topological sector is the construction of a 
height model on a bipartite lattice. Given the rules developed by Henley [4], a function 
h(x, y) can be constructed to assign an integer height to each vertex on the lattice dual 
to the one on which the dimers reside. The construction occurs as follows: set the height 
to be zero at an arbitrary plaquette, then assign integer heights by turning clockwise 
around the sites of the A sub lattice, adding +1 when crossing a bond with no dimer 
and subtracting z — 1 when crossing a dimer. Here z is the coordination number of 
the dimer occupied lattice. For the honeycomb lattice the dual lattice is the triangular 
lattice and assigning heights to its vertices amounts to assigning integer values to each 
plaquette of the honeycomb as shown in Figure 4. 



Classical Topological Order in Kagome Ice 



7 



If the height mapping is constructed for both the dimerization of interest and 
for the reference state as in Figure 4 then one can examine the difference in heights 
8h = h\c) — h\c) mapped on the transition graph, for which the following holds: 8h 
is constant inside each loop of the transition graph, Sh changes by +z when crossing 
a clockwise loop of the transition graph, and the winding numbers correspond to the 
average height difference between both sides of the sample. This last statement allows 
recovery of the winding numbers via the height model by taking the difference of heights 
at the edges of the lattice and dividing by the coordination number (Equations 7-8). 
For a lattice with periodic boundary conditions the height difference is taken between 
what the next height on the lattice should be and what it is upon wrapping around the 
boundary: 

Sh 

Sh x = h(2L,y) - h(0,y), w x = — -, (7) 

z 

Sh 

6h y = h{x,2L)-h(x,0), Wy = ^. (8) 

Thus, all three methods of characterizing the topology of the groundstate (transition 
graphs, cut numbers, and height mappings) are related. In a numerical simulation, it is 
therefore only necessary to measure one in order to define the topological sector. In the 
rest of this paper, we employ cut numbers measured directly on the kagome lattice, since 
this avoids the need for constructing the dimer model on the honeycomb lattice. In the 
next section, we compare measurements of the topological sectors defined via the cut 
number to measurements of susceptibility using a Monte Carlo simulation procedure. 



3. Classical Monte Carlo and Measurements of Susceptibility 

We perform Monte Carlo simulations on the pseudo-spin kagome Ising Hamiltonian 
Equation (1). We used both a loop update algorithm [18, 19] (adapted to the problem at 
hand) and a standard Metropolis algorithm allowing single spin flips. The loop update 
begins by first choosing an "up" (a, = 1) pseudospin at random; the loop "head" 
then enters a triangle associated with this pseudospin, and exits through a "down" 
(<jj = —1) pseudospin. The algorithm is adapted to work both in the constrained 
(Wills) and unconstrained [111] cases. In both situations the algorithm searches for a 
pathway flipping up and down pseudospins sequentially. When a minority spin is to 
be flipped, the algorithm has no choice and progresses deterministically to the next 
triangle. For a majority spin, it chooses one of the two directions at random. For the 
unconstrained case, majority and minority spins change at random along the loop, while 
for the constrained case they are fixed everywhere by the broken Z2 symmetry. Since, for 
both cases the ground state consists of an extensive manifold of degenerate states, the 
loop algorithm works at T = where there is no internal energy, and where the specific 
heat goes to zero allowing us to assure ergodicity, or partial ergodicity among the low 
energy states. At T > 0, excitations appear above the degenerate manifold, leading 
to triangles which violate the local ice rules. If the loop encounters such a triangle it 
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aborts. Metropolis updates in conjunction with the loop moves ensures a stochastic 
evolution of the system in this case. 

We further allow two variations of the loop building algorithm. In the first we 
allow the loops to grow to any size including loops that span the periodic boundaries. 
We refer to this as the "long" loop algorithm. As demonstrated below, the loops that 
wind around the periodic boundaries of the system allow fluctuations of the topological 
sector. We also define a "short" loop, which is built using the same algorithm, however 
if the loop size exceeds some pre-determined number of sites Nfj^, the loop is aborted 
and the update is not performed. Detailed balance is not affected if this procedure is 
followed. Typically, we chose N™* < L where L is the linear size of the lattice (of 
N = L x L x 3 spins), in order to study the effects of short loops on the ergodicity 
between topological sectors. 

In addition to the usual thermodynamic estimators (such as energy and specific 
heat), we directly measure the topological cut number defined in the previous section 
(Figure 3). For simplicity, the cut numbers are referenced against a theoretical state \C) 
in which Hx, Hy = 0. Although this reference state is unattainable by construction, it is 
convenient in that it makes the cut numbers directly equivalent to the winding numbers 
w x ,w y . Our definition of the cuts means they are bound between < Hx,Hy < L 
rather than the original range of [-L,L] given for the winding numbers of Rokshar and 
Kivelson, as discussed more below. 

We now present simulation results for the topological cut numbers on finite- 
temperature Monte Carlo simulations of both [111] ice (h = 2 J) and Wills ice (h = 0). 
Finite-temperature annealing simulations were performed on Wills ice (which does not 
enter a topologically ordered phase) and on [111] ice. In the latter case, simulations 
were performed with long loops, and also with short loops restricted to iVj^f < L size. 

In Figure 5 we illustrate the average cut number as a function of temperature for 
Wills ice and [111] ice with short and long loops. In Wills ice the average cut numbers 
remain constant over the entire temperature range at a value consistent with a random 
and unconstrained distribution, that is with (H x )/L = (H y )/L = 1/2 . This null result 
clearly illustrates the absence of topological information in the unconstrained system, 
with the mean cut number remaining unchanged as one passes from the high to low 
temperature regime. For the constrained system the cut number evolves continuously 
on a temperature scale between 1 < T/J < 10 as the topological constraints are 
imposed. For the long loop simulations the mean cut numbers divided by L approach 
(H x )/L = (H y )/L = 1/3 and their values fluctuate around this mean. That is, the 
topological sectors fluctuate around a mean value consistent with broken Z2 symmetry, 
in which exactly 1/3 of the pseudo spins point up and the system is topologically 
constrained (but not topologically ordered). For the short loop simulation the behaviour 
is quite different: below the temperature scale on which the constraints are imposed the 
cut numbers freeze into fixed values, H x = 4 and H y = 5 in this case, and the system 
becomes topologically ordered. Restricting to short loops therefore leads to loss of 
ergodicity and topological order as the constraints are imposed. 
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Figure 5. Topological cut numbers on an L = 12 system for: (a) Wills Ice (h = 0), 
(b) [111] Ice (h = 2) with Short Loops, (c) [111] Ice (h = 2) with Long Loops. 



The cut number described above is a non-local quantity, measured directly by 
circumscribing the toroidal geometry of the simulation cell. By definition, topological 
sectors require non-local measurements in order to be quantified, but such measures are 
notoriously inaccessible to experiment. Hence one must address the question of whether 
local measurements can be indicators of topological order in this system. For this reason 
we now examine susceptibilities and their capacity to give information concerning the 
topological order of a constrained system. We first define a pseudospin magnetization: 



(m) 



and an associated pseudospin susceptibility, 

(m) 2 ]. 



N W 2\ 



(9) 



(10) 



Although, in this case, xp could be accessed through polarized neutron scattering on spin 
ice materials in the [111] geometry [15, 20] a more experimentally accessible quantity is 
the real spin susceptibility. The magnetization and susceptibility associated with these 
degrees of freedom are defined: 



(M a ) 



where 



a = x,y,z\ 
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Figure 6. (a) Pseudospin susceptibility, and (b) Real spin susceptibility. Each case is 
for an L — 12 system. Results for unconstrained Wills Ice are in black, [111] long loop 
spin ice in red, and [111] short loop spin ice in blue. 



Xr = -[(M-M)-(M} 2 ]. (12) 

Note the Xr can equivalently be written directly in terms of the microscopic degrees of 
freedom: 

Xr = ^E[< S *' S i>-< S *>< S i>] ( 13 ) 

i,j 

We simulate both the pseudo and real spin susceptibilities, for both Wills and [111] 
ice, using the same procedures as for the cut numbers, described above. Pseudospin 
susceptibilities are illustrated in Figure 6(a). Here, one sees that T\ P reaches a 
limiting value for both Wills and [111] kagome ice for T < 1, which corresponds to the 
temperature where the system settles into the disordered groundstate manifold. The 
low-temperature value of T\ P is zero for [111] kagome ice, indicating that fluctuations in 
pseudospin magnetization are frozen out. This result is consistent with the topological 
constraints being enforced everywhere on a temperature scale less than h, so that once 
the number of triangles not satisfying the uud constraint becomes exponentially small, 
the fluctuations in m fall exponentially to zero. This is true for both the long and short 
loop algorithm and so is a measure of the applied topological constraints, but is not a 
measure of topological order. This result should be contrasted with that for Wills ice 
whose pseudo spin susceptibility continues to fluctuate at low temperature, illustrating 
that it is exploring the full low temperature phase space of upp and uud states. 

Figure 6(b) illustrates T\ r as a function of temperature, which also shows different 
behaviour as T — > for our different models. Wills ice and the long loop simulation for 
[111] both approach finite values as T approaches zero. There is, however a most striking 
difference between long and short loop simulations for [111] ice. For short loops, T\r 
drops abruptly to zero for T « 0.5. From Figure 5 one can see that this drop corresponds 
to the onset of topological order. As the sectors can only change by introducing a loop 
that spans the periodic boundaries, when such loops are suppressed the sector becomes 
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Figure 7. Measure of the real susceptibility in [111] ice for system sizes from L = 32 
to L = 3 at a constant temperature of T = 0.01. The temperature multiplied by the 
real susceptibility is plotted against the inverse linear system size. 



frozen. This sector freezing manifests itself in the magnetic fluctuations, as short closed 
loops inside the system carry no magnetization and contribute zero to Tx r , so that it 
falls abruptly to zero as the concentration of topological defects disappears. System- 
spanning long loops do however carry magnetization, and the finite value for T\ r at low 
temperature is consistent with the fluctuating topological sector at low temperature. 
Restricting loop length therefore leads to a dynamical phase transition to a state with 
long range topological order and this evolution is apparent in the real spin susceptibility. 

We can estimate the limiting value T\r as T — ^ by approximating the kagome 
lattice by a Hushimi tree of three-fold coordinated corner sharing triangles (see the 
Appendix) and summing the combinatorial series corresponding to the ground state 
manifold of states [21]. Such tree calculations have proved to be extremely accurate 
for the nearest neighbour spin ice model [21] despite the fact that the closed nearest 
neighbour pathways appearing on the d-dimensional lattice are neglected. As in the 
case of nearest neighbour spin ice and the six vertex model on a square lattice [22], 
the tree calculation yields T\ T {L = oo) = 2, while finite-size scaling (Figure 7) at 
T = 0.1 for [111] ice with long loops indicates that T\r{L = oo) 2.2. For the 
constrained system, the closed pathways seem to have a considerable influence on the 
low temperature susceptibility. Interestingly, this is not the case for Wills ice: here, the 
tree calculation also yields the modified Curie law, Tx r (L = oo) = 2 and simulation 
yields the same result within numerical accuracy. This agreement is consistent with the 
absence of non-local correlations in the unconstrained system. 

Finally, we can examine the effect of restricted loop size on the real susceptibility 
of a subset of the topologically ordered system. Since we have seen that restricted local 
dynamics are a critical ingredient for the non-ergodicity between topological sectors, 
we can examine the effect of loop size on subset regions of varying size. In doing 
these simulations we are able to test the effect of different boundary conditions, as small 
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Figure 8. Measure of the real susceptibility in [111] ice for for a varying subset £ of a 
L = 32 lattice. This data is taken from several simulations which were identical except 
for their loop algorithm, which was limited to produce loops of a specific size in each 



windows in a simulation of much greater size, approach the situation of open boundaries. 
Performing simulations on a L = 32 system, we measure Xr on a subset of spins with 
size £ x £ x 3 at T = 0.01. From Figure 6 for loops restricted to any size smaller than 
L, Txr = and the system will be in a topologically ordered phase. In Figure 8 we 
show susceptibility results for loop size less than L, as a function of linear subregion size 
£. Interestingly, T\ r ^ for all £ < L, even if £ is larger than the maximum allowed 
loop size. Magnetic fluctuations are allowed within the widow of size £, even though the 
total topological sector is frozen, as short loops spanning the edge of the window carry 
magnetization. As £ gets larger than the maximum loop size, all magnetic fluctuations 
result from edge loops, hence T\r versus £ scales linearly to zero as £ — > L. As £ 
increases the results become independent of the maximum loop length, as long as this is 
less than L. This suggests that the major contribution to this measure comes from short 
loops, which is rather surprising given that the loop length distribution is expected to 
be power law up to the imposed cut off. In contrast, if the loop lengths are unrestricted 
the susceptibility of the window increases monotonically towards the fluctuating sector 
susceptibility as £ approaches L. These results illustrate the non-local nature of the 
topological information encoded into the system. The window susceptibility for the 
two cases is the same if £ <C L, but their relative difference becomes of order one as £ 
increases and the effect of the boundaries are felt over macroscopic length scales. 

4. Discussion 

In this paper we have employed Monte Carlo simulations using local and non-local loop 
updates to characterize the classical topological fluctuations occurring in the ground 
state of the constrained kagome ice model, realized by placing an external magnetic 
field along the [111] direction of the nearest neighbour spin ice model. We characterize 
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the topological sector in the groundstate using a dual-lattice dimer representation: 
topological numbers are defined by taking a "cut" across dimers directed from the 
A to B sublattice of the dual (honeycomb) lattice. Defined in this way, we show that 
Monte Carlo simulations using local dynamics (including size-restricted loop updates) 
loose ergodicity between topological sectors as the system settles into the groundstate 
manifold of states, leading to a topologically ordered state. Unrestricted loop updates, 
which allow loops to wind around the periodic boundary of the system, allow the system 
to fluctuate between topological sectors at all temperatures. 

Defining a susceptibility Xr using the 2D projection of the real spin ice moments, 
we show that the onset of topological order has signatures in Txr- Specifically, when the 
system is not allowed to fluctuate between topological sectors in the groundstate (i.e. 
in the presence of local updates), Txr = 0. If ergodicity between topological sectors 
is restored by allowing loop Monte Carlo moves of arbitrary length, Txr remains finite 
down to T — > 0. As the system settles into the ground state manifold, the susceptibility 
crosses over from Curie's law, x = 1/T to a modified Curie law reflecting the hidden 
non-local correlations, which we estimate numerically to be Tx r ~ 2.2. This is slightly 
above the value estimated from a calculation of triangular spin units on a Husami tree, 
Txr = 2, showing the importance of closed pathways of spins in constrained kagome ice. 

Throughout the paper we have compared and contrasted the physics of constrained 
([111]) kagome ice with that of unconstrained (Wills) kagome ice. We have illustrated 
how the constraints lead to the emergence of a hidden U(l) gauge field and to a highly 
correlated liquid state. If we take the full, unconstrained phase space of Wills ice, there 
is no emergent gauge field and the excitations out of the ground state manifold do not 
carry topological charge. In analogy with spin ice, monopole physics has previously been 
discussed in the context of artificial magnetic arrays with nano-magnets arranged on a 
honeycomb lattice [23, 24]. There, the phase space of orientations for the nano-elements 
is the same as that for the real spins in the kagome ice models discussed in this paper. 
However, in order to realize a gas of monopole quasi particles in two dimensions, one 
must have both an underlying emergent gauge field and effective Coulomb interactions 
between the topological excitations out of the constrained manifold of low energy states. 
The first criterion can therefore only be satisfied if the Z2 symmetry is broken and 
topological constraints are imposed. Although there is some evidence that dipole-dipole 
interactions between elements could spontaneously break this symmetry [25], it does 
not appear to be the case in experimental presented so far, with the result that the 
monopole description can only be valid in the low density limit, with excitations above 
a specific ordered array of nano-elements. 

In conclusion, we have shown that a precise definition of the topological number or 
sector must involve a non-local cut measure which circumscribes the toroidal geometry 
of the simulation cell. However, by measuring the susceptibility for spins in windows 
contained within a bigger system, we have demonstrated that Tx r can be used to 
determine whether or not a system is ergodic between topological sectors. By increasing 
the window size, or measuring spin fluctuations on increasing scales, on finds scaling 
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towards two well-defined limits, depending on whether topological sector fluctuations 
occur or not. Susceptibility measurements such as this may prove to be a valuable probe 
of topological order on classical systems such as spin ice in the future. 
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6. Appendix: Limiting Value of the Real Susceptibility 

Adapting the argument laid out in the thesis of Jaubert [21] for pyrochlore spin ice, it is 
possible to construct limiting values for the real susceptibility within the kagome plane. 
Assuming that the system approximates a grand canonical ensemble of infinite size such 
as a Hushimi tree, rather than a lattice of fixed size with periodic boundary conditions, 
one can proceed by expanding the sums in the expression for the real susceptibility as 
below. 



Xr — rp 



So) 



(14) 



Since (Si) = and (Sf) = 1 V i: 



Xr rp 



So) 



(15) 



The Curie Law Xr — ^/T, is recovered in the paramagnet phase, as when the 
constraints of the ice rules are broken there is no correlation between nearest neighbors. 
As increased thermal excitations cause a breakdown of the ice rules the Curie law thus 
dictates the limiting behavior for high temperature. To derive the limiting behavior for 
low temperatures we rewrite the sum in Equation 15 as an infinite sum over all nearest 
neighbors. 



E(s*-s ) = £ £ (s t -s ) 

n=l i£n th NN 



(16) 



Choosing an arbitrary spin in the plane to be So the dot product can be computed 
between So and its four nearest neighbors Sj. Treating the out of plane spin as fixed by 
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Figure 9. States that obey the ice rules within the kagome plane, with accompanying 
unit vectors di. 

the [111] external field, there are three potential spin configurations (see Figure 9) that 
satisfy the ice rules in the kagome plane, with the following expectation values. 



E<S,-So> = ^-^) = (17) 

T,(^-^) = k-\ + h = o (is) 

E(S,-S ) = i(i + i) = i (19) 

(20) 

The average value of all three configurations is: 

£(S, • S )all configs. = ~(0 + + 1) = ^ (21) 

Since there are four nearest neighbors for every spin the number of nth nearest 
neighbors is given by 4(2 n_1 ) = 2 n+1 . Taking the sum to infinity gives: 



OO OO 1 OO Q 

E^ = E2 n+1 (^r = Ei = 1 ( 22 ) 

n=l n=l n=l 6 

Leading to a limiting value of the real susceptibility identical to that found for the 
3-D Hushimi tree: 

Xr {T -> 0) ~ |. (23) 

This result also follows for Wills Ice, and can be derived through the same counting 
argument by ignoring the restrictions imposed by the ice rules. This leads to three 
more states of minimum energy, specifically the two- in/one-out states, contributing to 
the nearest neighbor inner product sum. However, the three new states have the exact 
same averaged sum as the spin ice accessible states. Thus the sum over all possible 
configurations and accompanying real susceptibility are the same as in [111] ice. 

In the Monte Carlo simulation, these results were found to accurately reflect 
the limiting value of the real susceptibility of Wills Ice, even for small system size. 
However, once the external magnetic field was added to the Hamiltonian the topological 
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constraints imposed on the boundaries via the ice rules prevented the system from 
reaching the theoretical T — > limit for [111] ice, regardless of system size. 



